Fluctuations of power injection in randomly driven granular gases. 
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We investigate the large deviation function 7Too(w) for the fluctuations of the power W(i) = wt, 
integrated over a time t, injected by a homogeneous random driving into a granular gas, in the 
infinite time limit. Our analytical study starts from a generalized Liouville equation and exploits 
a Molecular Chaos-like assumption. We obtain an equation for the generating function of the 
cumulants /i(A) which appears as a generalization of the inelastic Boltzmann equation and has a 
clear physical interpretation. Reasonable assumptions are used to obtain (J,(\) in a closed analytical 
form. A Legendre transform is sufficient to get the large deviation function 7100 (w). Our main 
result, apart from an estimate of all the cumulants of W(t) at large times t, is that 7r 00(10) has 
no negative branch. This immediately results in the failure of the Gallavotti-Cohen Fluctuation 
Relation (GCFR), that in previous studies had been suggested to be valid for injected power in 
driven granular gases. We also present numerical results, in order to discuss the finite time behavior 
of the fluctuations of W(t). We discover that their probability density function converges extremely 
slowly to its asymptotic scaling form: the third cumulant saturates after a characteristic time r larger 
than ~ 50 mean free times and the higher order cumulants evolve even slower. The asymptotic value 
is in good agreement with our theory. Remarkably, a numerical check of the GCFR is feasible only at 
small times (at most t/10), since negative events disappear at larger times. At such small times this 
check leads to the misleading conclusion that GCFR is satisfied for iToo(w). We offer an explanation 
for this remarkable apparent verification. In the inelastic Maxwell model, where a better statistics 
can be achieved, we are able to numerically observe the failure of GCFR. 

PACS numbers: 

I. INTRODUCTION 

In non-equilibrium statistical physics few general results have been established Q . A quick historical survey begins 
with the Einstein relation Q between the velocity fluctuations at equilibrium of a Brownian particle and the variation 
of its velocity under the effect of a small applied force. This relation between fluctuations and dissipation has been 
then explored and fruitfully extended in many ways, first by Onsager and then by Green and Kubo 0. All 
those works have been concluded with the establishment of the Green-Kubo formulae and the fluctuation-dissipation 
theorem, essentially a local formulation of the previous relations. Later, a thermodynamical phenomenological de- 
scription [|| was formulated, based on the hypothesis of local thermal equilibrium. Within the latter framework the 
entropy variations of a macroscopic system evolve according to the evolution equation 



<LS a- f dVV- J s • (1) 

J system 



dt 



The first term a is the irreversible entropy source, which vanishes at equilibrium, and is strictly positive out of thermal 
equilibrium. The second term J V • J5 is related to the external constraints (or fields), which are necessary to drive 
the system in a non-equilibrium state. In a steady state, of course those two terms are equal. Furthermore, as opposed 
to <7, whose origin is purely statistical, J ,5 is the average of a microscopically well defined quantity. We shall refer to 
this quantity as a mean entropy flux. In many physical systems this entropy flux is associated to a current, where the 
conveyed quantities arc usually the energy, the momentum, or the density. 

More recently, attempts to reveal the general behavior of large classes of systems have been performed, focusing on 
large deviation properties and on macroscopic quantities. For a long time the study of dissipative systems has been 
centered on local quantities, such as structure factors, correlations, and velocity fields. Nevertheless recent results Q 
support the intuition that in order to observe universal features, it is helpful to study "global" quantities. Averaging 
is indeed an effective way to bypass differences arising from microscopic details, keeping however essential features. 
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Most spectacular and surprising examples are provided in 0, , where it turns out that three different quantities 
(the total magnetization of the two-dimensional XY model near the critical point, the instantaneous injected power 
in a turbulent ffow, and the energy fluctuations of a cooling granular gas close to the clustering instability), suitably 
rescaled, display the same probability distribution function (pdf). 

In addition, we have learnt from equilibrium statistical mechanics that the most useful informations are contained 
in the state functions, the free energy being the main example. Those functions are, in the probabilistic interpretation, 
large deviation functions (ldf ) 1 . In equilibrium statistical physics, the large deviation principle enters when considering 
microscopic quantities summed over all the degrees of freedom, or integrated over the whole volume of the system. 
Thus, large deviations are associated to the size-extensivity of the macroscopic quantities of interest. As a matter of 
fact in non-equilibrium systems the dynamics, and hence the time, plays a central role (even in stationary states). 
This has led to the introduction of time extensive quantities, which can be easily obtained integrating a stationary 
quantity over a given time interval. The simplest choice to attempt for a universal description seems indeed to study 
the time large deviation function of a flux integrated over a time t. In this framework the infinite time limit could be 
seen as the "non-equilibrium counterpart" of the thermodynamic limit in equilibrium physics. 

With this objective Evans, Cohen and Morriss discovered a particular symmetry property of the ldf of the integrated 
injected power for a system of thermostatted sheared hard disks 9]. In order to have a non-equilibrium steady state, 
it is necessary to have a competition between at least two forces acting on the system. In the particular case of the 
energy injected by the shearing force is compensated by an external thermostatting force, created ad hoc, in such a way 
that the total kinetic energy of the system is constant. This symmetry property has subsequently been formalized in a 
theorem by Gallavotti and Cohen [l(| . The theorem has been proved for deterministic and time reversible dynamical 
systems, under the hypothesis of strong chaoticity. In this context the quantity S(t) is the phase space contraction 
rate integrated over a time t. It turns out that in a large class of dynamical systems the phase space contraction rate 
is endowed with the significance of the entropy flux defined in non-equilibrium thermodynamics. For thermostatted 
systems [ll| this is equivalent to the power injected by the thermostatting force divided by the temperature of the 
system The Gallavotti-Cohen Fluctuation Relation (GCFR from now on) concerns the large deviation function 

of the probability P(S, t) of a quantity playing the role of the entropy flux, defined as 

7Too(s) = lim 7r t (s), with 7r t (s) = ^ log P(st,t) . (2) 

t — >OD t 

Then, the GCFR reads: 

7Too(s) - 7Too(-s) = S . (3) 

Similar results were obtained for Markov processes by Kurchan ^| an( A by Lebowitz and Spohn [l3| . 

After the GCFR has been proposed, several attempts of verification in experiments and numerical simulations 
have been proposed, always with success |l lj . In many of these attempts however, the quantity S(i) was not the 
integrated phase space contraction rate, and the original hypothesis of the theorem were not guaranteed. The phase 
space contraction rate is indeed usually not accessible in experiments, and strong chaoticity (Anosov property) is 
not strictly satisfied in most real systems. The choice of the quantity S(t) usually goes to functional that are 
related to equilibrium or near-equilibrium entropy-like quantities. Typically one expects that the fluxes defined in 
non-equilibrium thermodynamics [fj are the macroscopic analogous of the phase space contraction rate, and this is 
definitely true in specific thermostatted models ^lj. The rate of heat transport from the thermostat, or equivalently 
the power injected by the thermostat into the system, divided by the temperature of the thermostat, is the natural 
candidate. In order to check for possible extensions of the GCFR for systems without time-reversal invariance the 
study of the pdf of the injected power seems the simplest choice. Such a program has been initiated by Aumaitre 
et al. 15], who measured the injected power distribution for several microscopically irreversible systems. Within 
the available accuracy, they concluded that the GCFR holds but importantly pointed out the possible relevance of 
spurious effects at large times where the necessarily limited data gathered may not allow to sample the regions where 
GCFR violations could exist. An analytical approach was also proposed by Farago 0, who analytically computed 
the ldf of a particular model without time-reversal symmetry, showing that GCFR does not hold, at least for this 
particular model. Nevertheless a recent experiment [17| observed the validity of the GCFR measuring the integrated 
injected power rendered dimensionless by an appropriate energy scale. 



1 As an example, for a system in the canonical ensemble with a free energy F{T, V, N) and a partition function Z = e @ F , all the physical 
features (such phase transitions, etc.) are contained in the free energy per particle / = F/N, which in the thermodynamic limit is 
/ = — A limjv^oo 1 °% Z i which only depends on intensive variables. 
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In a recent study, we have shown that in general the GCFR is not satisfied for the integrated injected power in 
homogeneously driven granular gases |l8l | and that the results of the experiment described in [T^ | were insufficient to 
claim its verification In this paper we discuss the derivation of an equation for the large deviation function for 
injected power in homogeneously driven granular gases, obtaining its solution under a series of reasonable approx- 
imations. The Gallavotti-Cohen Fluctuation Relation cannot be satisfied in this model because the injected power 
large deviation has no negative tail. We also provide details from numerical simulations of the model, showing how 
the GCFR may seem to be satisfied at small times [i.e. considering n t (s) instead of 7roo(s)] and small values of s. In 
experiments and simulations true large deviations are argued to be unreachable. 

In our approach, the granular gas is modeled as a dilute system of inelastic hard spheres. The heating mechanism 
that allows the gas to reach a steady-state is the widely used stochastic thermostat [2(J . It differs from the experimental 
setup of Ref. [l7| in that the energy injection is fully uniform (each particle is subjected to a random force) instead of 
being boundary driven. The model under study is however relevant for other two dimensional experimental situations, 
see e.g. |2^|. Our goal is to compute, analytically and numerically, the distribution of the power injected into the gas 
by those random kicks. 

In section ITT1 the inelastic Hard Spheres model is presented and the details of the derivation of the large deviation 
function of the injected power are presented. Section IlIII repeats the procedure for the closely related and much more 
tractable inelastic Maxwell Model (practical motivations will become clear later). Numerical results are provided in 
section llVl We then conclude the paper (Section 0) with a discussion of our results. 



II. THE INELASTIC HARD SPHERES MODEL 



A. The Model 



We consider a granular gas made up of N identical spherical particles undergoing dissipative collisions. The collisions 
between two particles i and j conserve the total momentum and reduce the normal component of the relative velocity, 
i.e. (v* — v*) -<r = — a(vj — Vj) -<r, where the stars denote the post-collisional velocities and <x the direction joining the 
centers of the colliding particles. Energy injection, achieved by means of random forces (kicks) acting independently 
on each particle, drives the gas into a non-equilibrium steady state. The equation of motion governing the dynamics 
of each particle is therefore: 

dt 1 1 ' W 

where F? oU is the force due to collisions and F* h is a Gaussian white noise (i.e. (F-^(t)F^g(t')) = 2r5ijSysS(t — t'), 
where the angular brackets denote statistical average; the subscripts i and j are used to refer to the particles, while 
7 and 5 denote the Euclidean components of the random force). This model is one of the most studied in granular 
gas theory and reproduces many qualitative features of real driven inelastic gases |20t |22| . After a few collisions per 
particle, the system attains a non-equilibrium stationary state. Furthermore, this state is homogeneous and does not 
develop spatial inhomogeneities, in contrast to what happens for example in the freely cooling state of a granular 
gas. From the equations of motion it is possible to derive the Boltzmann equation governing the evolution of the 
one-particle velocity distribution function / which, since the system is homogeneous, will not depend on the positions 
of the particles. If, in addition, we resort to the molecular chaos assumption, this Boltzmann equation reads 

9 t f(v 1 ,t) = J[/,/] + rA vl /( Vl ) , (5) 

where the Laplace operator A v = (d/d v ) 2 is a diffusion term in velocity space characterizing the effect of the random 
force, while J[f, /] is the collision integral, which takes into account the inelasticity of the collisions: 



Af, f] = \j dv 2 J' d<x(v 12 • «x) (^/(vD/(vr) ~ /(vx)/(v 2 ) j . (6) 

In the above expression of the collision integral the notation v 12 denotes the relative velocity between particles 1 
and 2, while the two stars superscript (i.e. v**) denotes the pre-collisional velocity of a particle having velocity v. 
Moreover the length I = (na d ^ 1 )^ 1 , where n is the particle density, a the diameter of the particles and d the space 
dimension, is proportional to the mean free path. Finally, the primed integral J will be used as a short-hand notation 
to denote J 0(vi2 • &), where is the Heaviside step function. The granular temperature of the system is defined as 
the mean kinetic energy per degree of freedom, T g = (v 2 )/d. Such a definition for the temperature is purely kinetic, 
and can a priori not be interpreted as a thermodynamical temperature. 
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B. The velocity distribution function 

In this section we briefly review the known results concerning the solution of the Boltzmann equation (jSJ). The 
stationary solution of this equation has extensively been investigated in the last years, but an exact solution is still 
missing. Nevertheless a general method is to look for solutions in the form of a Gaussian distribution multiplied by 
a series of Sonine polynomials: 

Mv)=e-&[l + f^a p S p (£r)) . (7) 



P =i 



The expression of the first three Sonine polynomials is: 

So(x) = 1 

S x {x) = -x+^d (8) 
S 2 {x) = \x 2 - i(d + 2)x + ^d(d + 2) . 
Moreover the coefficients a p are found to be proportional to the averaged polynomial of order p: 

a P = A p( S p{^)) ' ( 9 ) 

where A p is a constant. From this observation one directly obtains that the first coefficient a\ vanishes by definition 
of the temperature. A first approximation for the velocity pdf is therefore to truncate the expansion at second 
order (p = 2). An approximated expression for the coefficient ai has been found as a function of the restitution 
coefficient a and the dimension d [22l l23l l24j . It must be noted that this approximation is only valid for not too 
large velocities, since the tails of the pdf have been shown (2^ to be overpopulated with respect to the Gaussian 
distribution. It is known j2^| that at high energies log/(v) ~ ~(v/v c ) 3 ^ 2 with a threshold velocity v c that diverges 
when the dimension d goes to infinity. This means that at high dimensions the distribution is almost a Gaussian with 
small Sonine corrections. All the above results have been confirmed by numerical simulations, in particular through 
Molecular Dynamics (MD) and Direct Simulation Monte Carlo (DSMC) methods. Those two numerical methods, 
although very different, show a surprisingly good agreement. This points out the validity of the molecular chaos 
assumption and thus the relevance of the DSMC method, which is particularly well adapted to simulate the dynamics 
of a homogeneous dilute gas. 

C. Equation for the injected power 

In this section, we will show how to obtain a kinetic equation able to describe the behavior of the pdf of the 
integrated injected power at large times. The latter quantity is the total work W provided by the thermostat over a 
time interval [0, t]: 

W(t)= [ dt^Pf . Vi . (10) 
Jo 

Our interest goes to the distribution of W(t), denoted by P(W,t), and to its associated large deviation function 
■7Too(w) defined for the reduced variable w = W/t (W(i) being extensive in time): 

TtooH = lim 7r t (to) ,n(w) = -lnP(W = wt,t). (11) 

t—>co t 

We introduce p(Tn, VV, t), the probability that the system is in state Tn at time t with W(t) = W. The function we 
want to calculate is 



p{w.t) - / ar NP (r N ,w,t) . (12) 

We shall focus on the generating function of the phase space density 

p{T N , A, t) = [ dWe- xw p(T N , W, t) (13) 
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and on the large deviation function of 

P(X,t) = J dWc- xw P(W,t) = J dT N p(T N ,X,t) , (14) 



which wc define as 



p(X) = lim -lnP(A,i) . (15) 

i-^oo t 



Note that /x(A) is the generating function of the cumulants of W, namely 

(16) 



lim <^ = (-1)" 



i v ' dA™ 



A=0 



Moreover ir 00 ('w) can be obtained from p(X) by means of a Legendre transform, i.e. n 00 (w) = /u(A*) + A*w with A* 
such that //(A*) = — tu. 

The observable W is non-stationary and we are interested in the evolution equation for the extended phase space 
density p(T N ,W,t). It varies in time under the combined effect of the inelastic collisions (which do not trigger any 
change in W) and of the random kicks: 

dtp = d t p +d t p . (17) 



collisions 



kicks 



Instantaneous collisions do not modify the value of W, only the random kicks will affect it. In between two successive 
collisions, the amount of energy provided by the thermostat on each particle is AW = (v(t~ +1 ) 2 — v(i+) 2 )/2, where 
U denotes the time of the i— th collision and the — (+) superscript simply denotes that the value of the velocity at a 
collision time has to be taken before (after) the collision has taken place. 

In order to characterize the role of the thermostat on the generating function p it is useful to discretize the problem 
in velocity space, considering that the effect of the random kicks on one particle is a continuous time random walk 
of its velocity on a (I dimensional lattice of step a, with hopping rate 7. It then appears that the observable W is 
Markovian, and the master equation governing the evolution of the probability of being at site n (which represents 
the discretized velocity of one particle) at time t, with the value W(t) = n(t) 2 /2 reads: 

dtP(n, W, t) = 7 (P(n + a, W + n • a + a 2 /2, t) + P(n - a, W - n • a + a 2 /2, t)) - 2rf 7 P(n, W, t) . (18) 

a 

Here {a} is a set of d orthogonal vectors defining the lattice. This yields, in terms of the generating function 
P(n, A', t) = J dWc- Vw P(n, W, t), 

d t P(n, A', t) = jY ( p ( n + a > A '> t)e A ' n - a+A ' a2/2 + P(n - a, A', ()c- A '"' a+A ' a!/2 ) - 2d 7 P(n, A', t). (19) 

a 

Next, it is possible to expand the above expression in powers of a at order a 2 and consider the continuum limit a — > 
with the scaling variables v = no, T = ja 2 and A = X'/a fixed. Finally this yields, considering that the thermostat 
acts independently on each particle, 

dtp - V fr(A V( + 2AIV, • d Vi + T(dX + X 2 v 2 )} p (20) 

kicks — L J 

i 

This additional piece is linear in p, just as the collision part. The large time behavior of p is thus governed by the 
largest eigenvalue p,(X) of the evolution operator of p. In the large time limit, we expect that 

/ 5(r Ar ,A,i)^C(A)e^ A ) t ( 5(r 7V ,A), (21) 

where C(A) = 1/ J dTiypiTN, A), and p(Tn,X) is the eigenfunction associated to p, which for convenience has been 
taken normalized to unity (in order that C(A) ^ the initial state must have a nonzero projection onto this eigen- 
function). We then introduce 

fW(v,X,t) = J dT N - 1P , (22) 
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where J dT^—i means an integration over N — 1 particle coordinates. This function obeys 

d t f (1) (v, a, t) = rA v /w + 2Ara v ■ v/« + r(A V - + J (23) 

with J = J dWe~ xw J the Laplace transform of the collision integral J in which the usual velocity distributions 
/^(vi,t) and /( 2 ) (vi, v 2 , i) are replaced by /^(v, W, t) and /^(vi, v 2 , W,t) involving W as well as the velocities. 
Quite unexpectedly the above equation has a clear physical interpretation: consider a many particle system where a 
noise of strength T and a viscous friction-like force F = — 2ATv act independently on each particle, and where the 
particles interact by inelastic collisions. Consider then that the particles annihilate/branch (depending on the sign of 
A) at constant rate dXT, and branch with a rate proportional to A 2 i» 2 r. Then, the equation governing the evolution 
of the one particle velocity distribution of such a system is exactly the equation (|23(l . where A is a parameter tuning 
the strength of the external fields. In spite of there being no a priori reason for that, p, as well as /W = JdlV-ip, 
can consequently be interpreted as probability density functions. 

Since p(T N ,X,t) ~ C(X)e^ t p{T N , A), the one and two-point functions /W(v,W,i) and /< 2 ' (v 1; v 2 , W, t) that 
enter the expression of J are expected to verify in A— space, at large times, 

/W(v 1 ,A ) t) = C(A)e^/ (1) (v 1 ,A) ) (24) 

and 

/( 2 )(v 1 ,v 2 ,A,i) = C(A)e^/ (2) (v 1 ,v 2 ,A), (25) 

where both /W = j dTN-iP and = J <1Tn^2P are normalized to unity. We perform the following molecular- 
chaos-like assumption: 

/( 2 )( Vl) v 2) A)~/«(vx,A)/W(v 2 ,A) (26) 

which does have a definite physical interpretation in the language of the inelastic hard-spheres with fictitious dynamics 
(viscous friction, velocity dependent branching/annihilation) described in the above paragraph. Then we get that 



/i/(v, A, t) =TA v f + 2AIV • d v f + T(dX + A V)/ 



dv 2 d(TVl2 • <x 

V12 -CT>0 



«- 2 /K* , A)/(vr, A) - /(v x , A)/(v 2 , A) 



(27) 



where we have now omitted the superscript (1) denoting the one-point function. Note that the truncation 1261) 
respects two physical requirements. First, it does conserve W throughout a collision. This can be observed rewriting 
the product of the two one-point functions as 

/ (vi** } , A) / (v 2 **> , A) = J dW dWi dW 2 5(Wi +W 2 - W)c- xw g (v^ , Wi) g (v 2 " } , W 2 ) , (28) 

where g(v,W) — J dr^y-i p(Tn, W). This leads to the simple remark that if the particle 1 carries a fraction Wi 
of the integrated injected power before a collision, and the particle 2 carries a fraction >V 2 , then the global amount 
W = Wi + W 2 is preserved in a collision. Second, the A = limiting case yields the usual Boltzmann equation, since 
in this case a stationary solution exists, and hence p(X = 0) = 0. The boundary condition to the evolution equation 
above is thus: 

/(v,A = 0) = / st (v) (29) 
with /st(v) the stationary velocity pdf (cf. Eq.lQ)). 

D. The limits for large velocity and for A — > oo. 
Equation (|27|l can be written in the form: 

A«/(vi , A) = Te- A V1 ( /( Vl ,X)]+J(f, /) . (30) 
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Defining a new function F(y, A) = exp f^f-^/(v, A) the above equation reads: 

/i(A)F(v 1) A)=A VlJ F(v 1 ,A) + J X [F,F] , (31) 



where 



Jx[F,F} = J dv 2 J' d&(v 12 -&)e~^ ^e-^(^- 1 )( Vl2 ^ 2 F(vr,A)^(vr,A)-F(v 1 ,A)^(v 2 ,A)j (32) 

is the collision integral. Repeating the arguments by van Noije and Ernst we shall now compute the behavior of 
the tails of /. In the large velocities limit one has v±2 ~ v%, and for A > the gain term of the collision integral can 
be neglected, since the argument of its exponential prefactor is strictly negative for inelastic collisions (a < 1). In 
this limit the solution of eq. (|31|l is of the form: 

1/2N 




log/(v, A) = -A- - ~-\ ^ - ^ fcvV* + oV) , (33) 



where (3\ — ^ d 1 ^ 2 /r((d + l)/2) comes from an angular integration. This result is important since it shows that 
the high velocity tails of the function /(v, A) are Gaussian. For A = the known result concerning the non-Gaussian 
tails of the one-particle velocity pdf |23 is recovered. 

Moreover, in the limit A — > oo the gain term of the collision integral can be neglected. Thus, the resulting equation 
can be interpreted as a Boltzmann equation for a system of particles subjected to a random force, interacting with 
an "effective potential" which makes the collision kernel oc (vi2 • <x) exp — Airf/2, and annihilating when a collision 
occurs. The Boltzmann equation of this system reads: 



8 t F = A V F - J dv 2 / dd-(vi2 • &)e~r F( Vl )F(v 2 ) . (34) 

This immediately shows that /x(A) for high values of A is always negative, since the density of particles (given by the 
integral of F) is clearly decreasing. The consequence of this observation, together with the fact that fj,(X) must be 
always convex and is zero in A = 0, is that fjf(X) < for any value of A . For negative w, there is therefore no possible 
A* satisfying //(A*) = — w: the immediate consequence is that ^^(w), i.e. the Legendre transform of /x(A), is not 
defined for negative w. As announced in the introduction, this is a key point in invalidating the GCFR. 



E. Energy balance 



Before going further in the analytical calculation of ttoo(u>), we provide a complementary argument for the absence 
of negative tail. The reasoning takes as starting point the fact that the energy variations in a time t are given by the 
difference of the injected power integrated over a time t, W(i), and the energy dissipated through collisions over the 
same time interval T>(t). Hence the fluctuating total kinetic energy E(t) = J^. vf/2 varies according to 

AE(t) = E{t) - E(0) = W(t) - V(t) , (35) 

where T> > 0. It is important to note that the above equation is very general, and applies to many dissipative systems 
where an external driving supplies the energy lost by some dissipation mechanism in order to reach a non-equilibrium 
steady state. Let us define the probability V(z) of having AE(t) = z. If the total work W(i) is sampled starting from 
a fixed initial condition E(0) = const, then it is easy to show that the left tail of the distribution V is bounded by 
— E(0) (since E(t) > 0). This implies that the large deviation function associated to the stationary quantity AE(t) 
has no negative contributions, as well as the large deviation function tt^ of the total work W = T> + AE, since it is 
the sum of two quantities with no large negative events in the long times limit. On the contrary, if the total work 
distribution is obtained sampling over segments of a unique trajectory, the above argument does not hold anymore, 
since both E{t) and E(0) are fluctuating quantities. Thus, the pdf V(z) is no more bounded, but, for times much 
longer than the characteristic energy correlation time, it is symmetric with respect to the z = axis. Our purpose 
in this section is to give some phenomenological arguments to understand under which conditions on V(z) it is likely 
that W and T> have the same large deviations function. Since AE is not a time-extensive variable, in principle it 
does not have large deviations. Nevertheless a large fluctuation of AE may affect the behavior of the time-extensive 
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quantities W or V. Hence we will focus our attention to the tails of AE, which we will suppose to have a stretched 
exponential behavior with an exponent S (i.e. V(z) ~ exp — \z\ s when z — > ±00). It is useful then to note that 

( f , if S < 1 

C(e) = lim -log7>(A.E = et) = I -lei , if S = 1 (36) 

t— >oo t 

y 00 , if 5 > 1 (and e ^ 0). 

We therefore observe that the exponential distribution (5 — 1) plays a limiting role for the function £. Besides, if 
one looks at what happens in Laplace space (i.e. to the large deviation function of the Laplace transform of V\ the 
infinite time limit does not depend anymore on the particular value of the exponent S: 

0(A) = lim - log(c AA£; ) = . (37) 

Now a natural question arise: is it possible to recover the function £ knowing ? Usually the Legendre transform 
is the link between those two functions. It is easy to see that the Legendre transform of = gives £ = —00, 
independently of 5. As we will show, this "failure" of the Legendre transform is related with analyticity breaking of 
the function 0. If S > 1 the Laplace transform of V is analytical in the whole complex plane of A. On the contrary, 
if S = 1, the Laplace transform of V is very likely to presents singularities, as well as analytical cuts in the complex 
plane of A. The Legendre transform as an inversion formula of the Laplace transform is obtained using the saddle 
point method when integrating the inverse Laplace integral on a straight line parallel to the imaginary axis of the 
complex A plane. It is seen indeed that if (exp XAE) presents singularities or cuts, then the Laplace transform cannot 
be inverted anymore. 

Returning to our original problem, it now seems reasonable to argue that in a general way the large deviations 
functions associated to the Laplace transform of the pdf of W and T> are the same (i.e. lim(exp(— XW))/t = 
lim(exp(— XD))/t). However, if T'(AE) has exponential tails, non-analyticities may arise in Laplace space, and the 
large deviations functions of W and T> could be different. Besides, if V(AE) has tails decreasing faster then the 
exponential, the large deviations functions of W and V are the same. 

This scenario gives then simple hints for the knowledge of the large deviation functions of two time-extensive 
quantities which differs only by a stationary term. It turns out that if this stationary term is distributed with 
exponential tails, its fluctuations may affect the large deviation functions of the time-extensive quantities. This is 
clearly already seen in eq. I|36|l . which shows that in this case fluctuations of order t are possible. This clearly does 
not happen if the tails of the pdf of the stationary term decrease faster than any exponential. In this case fluctuations 
of order t vanish in the infinite time limit, and the two large deviation functions are equal. Moreover, this scenario is 
perfectly in line with the explicit results obtained by Farago in several simple models |0, |2(| , and by Van Zon and 
Cohen in |27|. In all those solvable cases the system considered is a one-particle system, with a very low number of 
degrees of freedom. Hence, if the velocity of the particle has a Gaussian distribution, the energy pdf has exponential 
tails, and in principle the distribution of W and V can be different. Nevertheless, if the system under investigation is 
a many particle system, as it is the case here, the energy pdf is distributed following the gamma distribution : 

P(E) ~ E% ~ l exp(-E) (38) 

where N is the number of degrees of freedom of the system. Moreover, in the thermodynamic limit (i.e. when ./V — > 00) 
the above expression of P(E) approaches the Gaussian distribution. Thus, taking the thermodynamic limit before 
that of large times, it follows from the above discussion that for large systems (i.e. in the thermodynamic limit) the 
total work VV and the dissipated energy T> have the same large deviation function. This is in line with the observations 
of the previous subsection, which qualitatively predict the absence of a negative tail for P(W). 



F. The cumulants 



In this section we find an approximated expression of /i(A) solving a system of equations obtained projecting Ij27(l 
on the first velocity moments. First we shall define a dimensionless velocity c = v/i>o(A), where t>o(A) plays the role 
of a thermal velocity: 

v 2 (X) = 2T(A) = ljdw 2 /(v, A) . (39) 
Then, defining the function /(c, A) = vo(\) d /(v, A), and its related moments of order n 

m„(A) = / dcc"/(c,A) , (40) 
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one obtains the following recursion relation: 

r 



((J, + T(2n + d)X)m n — — n(n + d — 2)m„_ 2 + TX v rn n+2 — v$v n , (41) 



where 



u n = -Jdcc n J[f,f]. (42) 

Recalling the definition of the cumulants H16(l . and the approximated solution for the stationary velocity pdf, it 
appears natural to argue that, for A ~ 0, the function /(c, A) should be well approximated by an expansion around 
the Gaussian, in Sonine polynomials: 

/(c, A) - 0(c) (1 + a 1 (X)S 1 (c 2 ) + a 2 (X)S 2 (c 2 )) + 0(a 3 ) , (43) 

where (f>(c) = n~ d / 2 exp(— c 2 ) is the Gaussian distribution. Even in this case, from the relation © and from the 
definition (|39[) . the coefficient a\ is found to be 0. The method consists in taking the equation l|41|l for n = 0, 2 and 
4 in order to find an explicit expression of fi, vq, and a 2 in the limit A — > 0. The quantities v 2 and v± have been 
calculated at the first order in a 2 [22^ . and their explicit expressions are: 

(1-a 2 ) n d f 3 1 dT f 3 1 
^2 = V og 1 JL {l + -a 2 \ = -== 1 + —a 2 , (44) 



and 



with 



2< 16 Z J v /27f I 16 



^=^={T 1 + a 2 T 2 }, (45) 



2i = d + I + a 2 (46) 



2 



T 2 = — (10 d + 39+10 a z ) + j- (47) 

32 (1 — a) 

2/3 

where To = ( ttz^ttt ] is the the granular temperature obtained averaging over Gaussian velocity pdfs (i.e. the 



zero-th order of Sonine expansion). The expression of the first moments m n is: 

rn = 1 (48a) 
to 2 = d/2 (48b) 

to 4 = ( 1 + Q2 ^ 2 + d > (48c) 

to 6 = (1 + 3**) <*(2 + <0 (4 + d) 



With the help of the above defined temperature scale To, we shall now introduce some dimensionless variables: 

(49) 



M = A 1 ^ ? A = AT , 



v 2 - v — ^EKv 

Note that this scaling naturally defines the scales for the other quantities of interest, namely: 

T w ~ W 

* = W= dT> W = W)- m 

The expression of the moment equation (|41|l becomes, for the above defined dimensionless quantities: 

(jid + (2n + d)X^j m n = — ^ — — to„_ 2 + 2X 2 v 2 m n+2 - v i> n . (51) 
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First we solve the above equation for n = 0, getting the following result: 

A(A) - -A + A 2 v 2 (A). (52) 

Recalling that when A — > one has v 2 = 2T g + 0(A), it is important to note that if we restrict our analysis to the 
Gaussian approximation for P(W,t), that is if we truncate /x(A) to order A 2 , Eq. I)52|l will read: 



Then we see that indeed 



(/r - A(AT S - 1) . (53) 



/i(A)=/i(^-A) , (.vl :) 



T 9 



which means that -k^w) = max>{/i(A) 4- Xw} verifies 



w 

^oo{w) - TTcvi-w) = — . (55) 

This is, up to a prefactor T g , the GCFR. However, the nontrivial functions m„(A) will break the property l|54|) . as we 
shall explicitly show later. In order to characterize more precisely the dependence of jx upon A for small values of A , 
it is useful to expand Vq and 02 in powers of A: 

v 2 (\) = vf } + Xvf + X 2 v 2<2) + 0(X 3 ) (56a) 
o a (A) = 4 0) + Xa { 2 1] + X 2 a {2) + 0(A 3 ) (56b) 

In this way we can find Vq ( ' ( a 2*^J solving equation l|51|l for n = 2: 



(57) 



vf^-U^-^, (58) 



^ 2, =2-4 0) (i + ") + ^-^- (59) 



2 _ a 2 

12 ' 3 J ' 3 8 

Then we substitute Vq(X) in the third equation and expand it in powers of A to find the expression oiaf{a). Note 
that one has also to expand in powers of a-i and keep only the linear terms in order to be coherent with the v p -s 
calculations. We find the following expressions, which are plotted in Fig. 2] 



,(0) 



4 (1-a) (1-2 a 2 ) 



19 + 14 d - 3 a (9 + 2d) + 6(1 - a) a 2 



(60) 



(i)^ 4 (1-a) 2 (-l + 2a 2 ) (31 + 2a 2 + 16d) 
° 2 (19+14d-3a(9 + 2d) + 6(l-a)a 2 ) 2 1 ' 



(2) _ A(a) 
B(a) 



= (62) 



with 



A(a) = 16 (-1 + af (-1 + 2 a 2 ) {906 + a [-984 + a (85 + 3 a (-19 + 6 (-1 + a) a))] + 985 d+ 
+ a [-951 + a (-25 + 3 a (7 + 6 (-1 + a) a))} d + (269 + 3 a (-75 + 2 a (-7 + 3 a))) d 2 } , 




FIG. 1: and versus a for inelastic hard 

spheres driven by random forces in d — 2. 



FIG. 2: The solid line shows jl in the limit d — * oo for 
inelastic hard spheres. The dashed line is jl at fourth order 
in A from I52H for d — 2 and a = 0.5. Finally the dotted 
line shows the same quantity calculated with a truncation 
at second order in A, which would satisfy the GCFR. 



and 

B(a) = 3 (-19 - 14 d + 3 a (9 + 2 (-1 + a) a + 2 d) f (64) 

The_Un <0) expression, as well as the a 2 ' expression, coincide with the usual results established for granular gases 
|22ll23| . At this point the computation of the cumulants becomes straightforward. From relation l|16fl it follows: 

lim = (-l)»drT?- l n\ if"' . (65) 

Moreover, since the a!j corrections are numerically small, as shown in Fig. ^ the zero-th order (Gaussian) approx- 
imation already gives a good estimate for the cumulants. Namely, the first cumulants are , in this approximation: 



(W) c = tNdT , (W 2 ) c = 2tNdTT , 

(W 3 ) c = StNdTT^ , ( w4 ) c = -18tNdTT3 . ^ 

All the above expansions in powers of A, at the second order in Soninc coefficients (e.g. 02) can be carried out by 
expanding vq and 02 in H56|) to higher powers of A. Moreover, expanding in higher order in Sonine coefficient (e.g. 03) 
remains in principle still possible, but will involve a higher number of equations in the hierarchy l|51|) (e.g. n = 6), 
and therefore will need the expression of higher order collisional moments (e.g. vq). 



G. The solvable infinite dimension limit 

The previous discussion strongly suggests that at high dimensions /(v, A) is not far from a Gaussian. First, the 
inspection of the d — > 00 limit of the Sonine expansion performed in the previous subsection indicates that the 
corrections to the Gaussian (see formula IrjUl and loT)l carry a relative d _1 factor, thus leaving the Gaussian pdf as the 
leading behavior when d — > 00. Second, the A = case is well known to reproduce a Gaussian in the d ~ * 00 limit: 
there is numerical evidence that the collision integral is solved (i.e. is set to zero) by a Gaussian /, when d — > 00. Of 
course if the Gaussian solves the collision integral then it solves the whole equation. 

We are therefore led to consider, in the limit d — > 00, /(v, A) to be a Gaussian with a A-dependent second moment. 
In this situation the dimensionless function / will read: 



(67) 
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with c(A) = v/wo(A). In this context one can solve equation i|51|) in order to get an explicit expression for /i(A). The 
equation defined by (|51l) for n = gives: 



A(A) = -A + A 2 £ 2 (A), 
where Vq(X) is obtained from Eq. I|51|l for n = 2, which reads: 

* 2 «o - «o - 2A ^o + 1 = 0. 
The unique solution of the above equation which verifies the physical requirement Vq(0) = 1 is: 

2 



(68) 



(69) 



1 + 4 A 3 6i(A) 1 
4A 4 + ^~~2 



^ (1 + 4^ MA) 
A 2 2 A 8 y ' K ' 46i(A) 



(70) 



with 



& X (A) 




6 2 (A) + 6 3 (A), 62(A) 



A 3 (9 + \/3\/27 + 256 A 3 



(71) 



63(A) = 



9 + V3 V27 + 256A 3 



2\ 3§ A 4 



64(A) 



32 
A 3 



24 ( 1 + 4 A 3 



1 + 4 A : 



A fi 



A 12 



This expression of the velocity scale reduces to the kinetic temperature for A = 0, and decreases monotonically as 
A" 1 / 2 when A — > 00. This means that in the limit A — > 00 / approaches a Dirac distribution as exp(— Aw 2 /2). This 
feature supports the intuition that the small W events (which are related to the large values of A) are provided by 
the small velocities. The behavior of jl is shown in Fig. [21 The large deviations function p,(X) becomes complex 
for A < — because of the terms containing V27 + 256 A 3 . Moreover for large A the behavior of this function is 
/i(A) ~ —A3. In the vicinity of the singularity (i.e. A = Aq = — ^irs) the behavior of the large deviation function is: 



M(A) 



3 
2 2 7 3 



3 3/2 2 l/6^_ Aq + _ Aq) 



(72) 



From the behavior for large A it is possible to recover the left tail of the large deviation function Tr^. In general, 
if /i(A) ~ — \P for A — > 00, this leads to //(A») = — /?A* _1 = —w. This last relation tells us that for (3 < 1 we are 

recovering the limit w — > + , with a behavior of the large deviation function given by Troo(w) = /i(A») + \*w ~ wf- 1 . 
Moreover, from the behavior of \x near Ao, an analogous calculation provides the right tail of the large deviation 
function: "Koq(w) ~ Aqw, when w — > 00. Finally, in our particular case, the tails are given by 



TToo(w — > + ) — -W 1/3 , 7foo(u5 — > 00) — -W , 



Note that, as expected from the discussion in subsection III El there is no w < tail to tt c 
function tTqc^w) is depicted in Fig. [5] 



(73) 

The graph of the whole 



III. INELASTIC MAXWELL GAS 



A. The model 

The Maxwell gas is a kinetic model due to J.C. Maxwell, who observed that a pair potential proportional to r - 2 ( d - 1 ) j 
being r the distance between two interacting particles, gives rise to a great simplification of the collision integral |29| . 
In fact this kind of interaction makes the collision frequency velocity independent. It must be noted that when the 
inelasticity of the particles is considered, this model looses its well-defined physical interpretation, but it nevertheless 
keeps its own interest. The collision integral is analytically simpler than the hard particles model and preserves the 
essential physical ingredients in order to have qualitatively the same phenomenology. In the recent development of 
granular gases this kinetic model has been extensively investigated |30l l3ll l32l l33| . In particular Inelastic Maxwell 
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IB 




-0.5 



FIG. 3: Large deviation function 7Too(w) for d — oo, as obtained by numerical inverse Legendre transform of Eq. I|68|l . 

molecules display strong non-Gaussian velocity pdfs (generally stronger than for the hard spheres), which can be even 
pathological in some cases. For example the velocity pdf of the cooling regime of such a gas has tails with algebraic 
decay. Nevertheless for the stochastically driven gas such a pathology does not exists, and the tails of the velocity 
pdf have an exponential decay. Moreover, as it will appear in the next section, with this model the DSMC simulation 
algorithm is highly optimized, and thus a much more precise observation of rare events can be achieved. 
The Boltzmann equation governing the single particle velocity distribution is: 



d t f{w 1 ,t) = J M [/,/]+ rA vl /( Vl 

where Jm is the collision integral, which slightly changes from (|rj|): 



/(vi)/(v 2 ) 



(74) 



(75) 



The solution of such equation have been extensively analyzed |3(i l3ll l32| , and it turns out that even in this case a 
Sonine expansion up to the second order provides a good estimate of the pdf for low velocities, while the tails of the 
pdf are exponential. 

In order to compute the large deviation function associated to the total work distribution, the same reasoning 
carried out in the previous part still holds. Thus, the analogous of equation (|27|l becomes, for Maxwell molecules: 



/i/(v, A, t) =TA v f + 2Arv • d v f + T(d\ + A V)/ 

+ u; J dv 2 J da [^/(vr.^/cvr,^-/^!,^/^^) 

where the same molecular-chaos-like assumption as in section [H] has been implemented. 

B. The cumulants 



(76) 



In this section we will provide the computation of the cumulants for the integrated power injected by the thermostat. 
The procedure is formally the same as in the previous section, for the hard sphere gas. The recursive equation involving 
the velocity moments is found from (|76[) . and reads: 



(^ + T(2n + d)X)m n = —n(n + d- 2)to„_ 2 + rA 2 VQm n+2 - Lu v n , 



(77) 



where m n is the n— th moment related to the dimensionless function / (cf. eq 140|)). The expression of the collisional 
moments is now slightly different. The two first non- vanishing moments are |34j : 
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and 



»4 = ^-{Tx + a 2 T 2 ) , (79) 



with 

Ti = 2(5 + 3a 2 + Ad) , (80) 



T 2 = 3a 2 - °(17 + 4d)-3(3 + 4d) 
1 — a 

Here again T = — ^z^p. ) denotes the zero-th order temperature. Let us introduce the following rescaled variables: 

M = M^rp * = XT o 

rr (82) 
"° =2^ " p = ~W Up 

The equation for the moment reads for the dimcnsionless variables: 

(dp. + (2n + d)X)m n = —^n(n + d — 2)m„_ 2 + 2A 2 5 2 m Jl+2 — du n . (83) 
2u 

Expanding again the solution in a Gaussian distribution multiplied by a series of Sonine polynomials up to the second 
order one obtains, from the first three equations of the hierarchy Q83[l. a closed set of equations involving p, Vq and a 2 . 
Those equations can be solved in the vicinity of A = and in the limit of small non-Gaussianity (|a 2 (A)| <C 1), when 
linearized with respect to the coefficient a 2 . In this framework, which should be valid in the limit of small values of 
A, the solution is: 



M (A) = ^4±iA _ a2 + (-! + >^T4A + 2A(-2-A + VrT4A)) + ^ ^ 
2 4 v 1 + 4 A 



with 



02(A) = ' 6(1 ' Q) (1 + Q) == . (85) 

(l + a) (7 + 3 (a -2) a-4d) + 8 (a - I) (2 + d) V1 + 4A 

Expanding this expression in powers of A, it is possible to recover all the coefficients a 2 of the expansion l|5tj|) . The 
first coefficients are plotted in Fig. These results allows us to directly compute the cumulants from the derivatives 
of p: 



(W n ) c = dTT"- 1 -M- 
V ' dA" 



(86) 

A=0 



It turns out that the corrections coming from the Sonine expansion are small with respect to the Gaussian order, 
which therefore already provides a good estimate. The values of the first cumulants are: 

(W) c = tNdT , (W 2 ) c = 2tNdTT , 

(W 3 ) c = \2tNdYT 2 , (W 4 ) c = 120tNdTT% . ^ 



C. The solvable infinite dimension limit 



In this section we will provide an exact solution for the generating function /z(A) in the infinite dimension limit. As 
already noted by Ben-Nairn and Krapivsky, in this limit the velocity pdf approaches the Gaussian distribution |3l| . 
Therefore, we will show that a Maxwellian with a variance depending on A is indeed a solution of equation (JBJ. 
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FIG. 4: a,2 '\ and a^ 2 ' versus a for the inelastic 

Maxwell model in d = 2. 
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FIG. 5: The solid line shows Jj, in the limit d — ♦ oo for the 
inelastic Maxwell model, given in Eq. ggj. The dashed 
line is p, at fourth order in A from 1841 for d = 2 and 
a = 0.5. Finally the dotted line shows the same quantity 
calculated with a truncation at second order in A, which 
would satisfy the GCFR. 



Because of the convolution structure of the collision integral, a great simplification occurs introducing the Fourier 
Transform of the velocity pdf: 



F(k,A)= J dve ik v /(v,A). 
Assuming that F is an isotropic function of the wave vector k, the equation (|76|l reads: 

liF(k) = -Tk 2 F{k) - \TdF{k) - 2Arfc^|^ - A 2 r^^ + u (F(^k)F(^k)) & - u F(k) (89) 



where 



1 — ol \ \ „ / 1 + a 



£ = l - 1 1 - ( — ) I e , v=[ Hr- ) > (90) 



and the brackets with a subscript 9 denote an angular average 



o 

where B is the beta function. Injecting a Gaussian solution 4>{k) = e ~ into the above equation yields the following 
relation: 

fx = -Tk 2 - XdT + 2Xk 2 T + \ 2 TT - \ 2 Tk 2 T 2 + iv (e~ k2{1 -f )ST ^ _ Wq _ ( 92 ) 
Then, expanding the exponential term in a power series and using that 

when d — > oo, one sees that the above relation holds for each k, with: 

yu(A) = -XdT + X 2 TT(X) (94) 
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FIG. 6: noo(w) for the Inelastic Maxwell Model in the limit d — > oo. 



and, 



T{X) _ 1 + ST+UT, (96) 

where To is the granular temperature. Hence, the explicit form of the generating function is, using the rescaled 
variables defined in l|82|) is: 



m(A) = 2 ' ( 96 ) 

This function has the same behavior as its analogous for the IHS gas, since it presents a cut in the negative real 
axis (for A < —1/4), and it decreases monotonically (as —A 1 / 2 ) when A — > oo. Furthermore this expression is, up to 
irrelevant constants, exactly the same as the one found in a simpler system, the Ornstein-Ulhenbeck process, which 
has been extensively analyzed by Farago in 0, [2(| . The analytical expression of the large deviation function for the 
total work is easily computed as the Legendre transform of the above expression of /i(A). It reads: 

(w - l) 2 

froo(w =~ , - Q(«>) • (97) 

4:1V 

This function is plotted in Fig. One can easily see that even for this model there are no negative events in the 
large time limit. The left tail of the large deviation function decrease to — oo as —l/4w, while the right tail has a 
linear decreasing behavior. 

IV. NUMERICAL RESULTS 

In this section the results of numerical simulations of the two models (inelastic hard spheres and inelastic Maxwell 
model) arc presented with particular attention to the Fluctuation Relation for the injected power. We emphasize 
that obviously only w t (w) can be accessed numerically whereas the GCFR bears on ir^iw). This requires a precise 
discussion of finite time effects, that play a crucial role here. The main requirement to verify the validity of the 
Fluctuation Relation is a clean observation of a negative tail in the pdf of the injected power. This poses a dramatic 
limit to the time t of integration of W(i). In numerical simulations, as well as in real experiments, at times larger than 
a few mean free times the negative tail disappears. On the other hand, at times of the order of 1-3 mean free times, 
the Fluctuation Relation appears to be correctly verified for the inelastic Hard Spheres model and slightly violated 
for the inelastic Maxwell model. The measure of the cumulants, anyway, gives a neat indication of the fact that the 
time of convergence of the large deviation function is at least 10 times larger and that the true asymptotic is well 
reproduced by the theory exposed in this article. 

The stationary state of a driven granular gas, modeled by equation (J5J, under the assumption of Molecular Chaos, 
is very well reproduced by a Direct Simulation Monte Carlo |2j|, [2^ . As a first check of reliability of the algorithm, 
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FIG. 7: Probability density function of the injected power, 
p(w, i) = tP(W(t) = wt, t) with t equal to 1 mean free 
time. In all three cases the value of the restitution co- 
efficient is a = 0.9. Other parameters are a) N = 100, 
T = 0.5; b) N=100, T = 12.5; c) N = 200, T = 0.5. The 
dashed line represents a Gaussian with the same first two 
cumulants. These distributions have been obtained with 
~ 1.5 x 10 9 independent values of W(t). 
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FIG. 8: Ratio of P(VV, and a Gaussian with the same 
first two moments, for the same parameters as in figure [7] 
a) corresponds to N = 100, T = 0.5, b) to N = 100, 
T = 12.5 and c) to N = 200, F = 0.5. The range between 
the vertical dotted lines is the useful one for the check of 
the Gallavotti-Cohen relation. It can be noted that the 
strongest deviations from the Gaussian behavior appear 
outside of this range. 



we have measured the granular temperature T g and the first non-zero Soninc coefficient a 2 = ((v 4 )/(u 2 ) 2 - 3)/3. The 
measured granular temperature is always in perfect agreement with the estimate. The measured a 2 coefficient is a 
highly fluctuating quantity and its average is in very good agreement with the theoretical estimate. 



A. Inelastic Hard spheres 



In Fig. [7| the probability density functions p{w, t) = tP{wt, t) (only for t equal to 1 mean free time) for three 
different choices of parameters iV, V (at fixed restitution coefficient a) is shown. The values of the first two cumulants 
of the distribution and their theoretical values are compared in tabled with very good agreement. In the same table 
we present also the measure of the third and fourth rescaled cumulants. The formulae for the first few cumulants are: 



(x 2 ) - (x? 



(x 3 ) -3(x 2 )(x) +2(x) 3 

(x 4 ) - 4(x 3 ) (x) + I2(x 2 } (x) 2 - 6(x} 4 - 3(x 2 ) 2 



(98a) 
(98b) 
(98c) 
(98d) 
(98e) 



N 


r 


<w(t)>/« 


(mf)c/t 


NTd 


2ATdT 9 




t K(t) 


100 


0.5 


100 


20835 


100 


21052 


0.200432 


0.355176 


100 


12.5 


2500 


13019125 


2500 


13157900 


0.201739 


0.361635 


200 


0.5 


199.9 


42009 


200 


42120 


0.141589 


0.175454 



TABLE I: First cumulants, together with the skewness a(t) — (W(t) 3 ) c / ((W(i) 2 ) c ) and the kurtosis excess n(t) = 

(W(t) 4 ) c / ((VV(i) 2 ) c ) 2 of the distribution of injected work P(W,i), measured with t equal to 1 mean free time for differ- 
ent choices of the parameters. 

The comparison with a Gaussian with same mean value and same variance shows that the pdf P(W, t) is not exactly 
a Gaussian. In particular there are deviations from the Gaussian form in the right (positive) tail. This is well seen in 
Fig- El K must be noted that the largest deviations in the right tail arise at values of W(t) larger than the minimum 
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W(t) available in the left tail, i.e. they have no influence in the following plot of Fig. [5] regarding the Gallavotti-Cohen 
symmetry. 

In Fig. El the finite time Gallavotti-Cohen relation 7Tt(w) — ir t (—w) = fi e ffW is tested for the same choices of the 
parameters. The relation, at this level of resolution and for this value of the time t (1 mean free time), is well satisfied. 
Moreover table ITU shows that the value of /3 e // is well approximated by 1/T g , as expected if the truncation of /x(A) at 
the second order were valid, see eq. (I55|) . In Fig. EUthe same relation is checked for different, slightly larger values 
of t (i.e. up to t equal to 3 mean free times). No relevant deviations are observed as t is increased. Moreover this 
figure is important to understand the dramatic consequences that larger t have on the possibility to check numerically 
the Gallavotti-Cohen symmetry: as t is increased, events with negative integrated power injection become rarer and 
rarer. This eventually leads to the vanishing of the left branch of P(W, t). 
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FIG. 9: Test of the Gallavotti-Cohen relation Tr t (w) — Tv t (—w) — f3 e ffw, for the same choices of the parameters as in Fig. 
The values of the slope /3 e ff of the fitting dashed lines are in table ITT1 



N 


r 




l/T a 


100 


0.5 


0.0100 


0.00955 


100 


12.5 


0.000402 


0.000382 


200 


0.5 


0.00995 


0.00952 



TABLE II: Factor of proportionality in the "Gallavotti-Cohen" relation compared with /3 — 1/T g 



The main conclusion is that no appreciable departure from the A 2 truncation is observed at this level of resolution. 
Much larger statistics are required to probe the very high energy tails of p(w,t). Further numerical insights make 
evident that the small times used to check the GCFR (t smaller or equal to 3 mean free times) are far from the time 
where the asymptotic large deviation scaling starts working. In Fig. 1111 indeed, we show the numerically obtained 
third cumulant of W(i) (rescaled by the first cumulant), as a function of time. The time of saturation is of the order 
of ~ 50 mean free times. The saturation value is in very good agreement with the value predicted by our theory 
cq. i|(jt)|) . Note that this value is not at all trivial, since the third cumulant for a Gaussian distribution is zero. For 
this value of t, the numerically measured TTt(w) is shown in Fig. ^1 rescaled by (w). The accessible range of values 
from a numerical simulation is dramatically poor and it is already remarkable to have obtained a good measure of 
the third cumulant with such a resolution. 

To understand the reason for a verification at small times of the GCFR we build upon an argument put forward 
in [l5j where it was noted that for small arguments, n can be linearized, resulting in a linear behavior of 7Tt(w)— 7tt(— w). 
To go further and compute the prefactor of the linear term, we remark that near w = the pdf of w is almost a 
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FIG. 10: Test of the Gallavotti-Cohen relation m(w) — 7Tt(— w) = j3 e ffW, for the system with N = 100 and T = 0.5 for different 
values of t. We recall that in this case (w) = 100. The dashed line has slope /3 = 1/T g . In the inset the corresponding p(w,t) 
are shown. 
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FIG. 11: Dimensionless third cumulant (W 3 ) c = 
(W 3 ) c /((W}Tj ) for several times of integration. The time 
is in units of the mean collision time. Note that the time 
at which a stationary value of the rescaled cumulant is 
reached is much larger than the characteristic time of the 
system (the collision time). The horizontal dashed line 
shows the analytical prediction of Eq. I|6fi[l : (VV 3 ) = 8. 



FIG. 12: Numerical measure of tt* for a time of 50 collisions 
per particle (when a stationary value for the rescaled third 
cumulant is reached), shown with circles. The full line 
corresponds to 7Too(if) encoded in Eq. (I68H and displayed 
in Fig. 01 



Gaussian. In the Gaussian case we immediately get 7T t (w) — 7r t (— w) — PeffW with f3 e ff = 2(W(i)) /(W(t) 2 ) c . The first 
two cumulants at small times are easily obtained considering an uncorrelated sequence of energy injection, obtaining 
(W(t))/t = NTd and (W{t) 2 ) c /t = (Q^Ff ■ v,) 2 ) c = 2NTdT g . The value of the slope f3 eff = 1/T g is a direct 
consequence of these simple relations. In this case the GCFR observed is the Green-Kubo-like (or Einstein) relation, 
which has been investigated numerically in |35| and analytically in |3(j , with the conclusion that such a relation holds 
within numerically accessible precision although it is, strictly speaking, invalid |36j |. Small deviations from a Gaussian 
appear, in first approximation, as small deviations from the slope 1/T g , but the straight line behavior is robust since 
the first non-linear term of nt{w) — ir t (—w) is not w 2 but w 3 . 
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B. Inelastic Maxwell Model 



Numerical simulations of the Inelastic Maxwell Model have been performed with a Direct Simulation Monte Carlo 
analogous to the one used in the Hard Spheres model. Thanks to the simplifications present in this model, we are able 
to improve the number of collected data by a factor of more than ten. The distributions of the injected power p(w,t) 
are shown in Fig. 1131 for some choices of the restitution coefficient a. The driving amplitude T has been changed in 
order to keep constant the stationary granular temperature T g . In Fig. 1141 we have displayed the deviations from the 
Gaussian of P(W, t). The non-Gaussianity of P(W, t) is highly pronounced, but again it is striking only in the positive 
branch of the pdf. We have tried, with success, a fit with a fourth order polynomial. We observe that the reliability of 
such a fit increases with a, which is reminiscent of Sonine corrections phenomenology for the single particle velocity 
distribution alluded to earlier. 

Finally, in Fig. El we have attempted a check of the Gallavotti-Cohen fluctuation relation. The relation seems to 
be systematically violated. This appears in two points: 1) the right-left ratio of the large deviation function is not 
a straight line; 2) the best fitting line has a slope which is larger than (3. The "curvature" (and the deviation from 
the slope (3) increases with decreasing values of a, indicating that the inelasticity is the cause of the deviation from 
the Gallavotti-Cohen relation. It should be noted that to achieve this result we have collected more than 4 x 10 10 
independent values of Wit), so that the statistics of the negative large deviations could be clearly displayed. 

We also note that this numerically accessible violation of the GCFR is of a different nature than the one analytically 
shown in the previous sections since it occurs at the finite times, i.e. when wt still displays negative w events. 
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FIG. 13: p(w,t) = tP(wt,t) for different values of a (at fixed constant temperature T g ) in the Driven Inelastic Maxwell Model 
measured at a time t equal to 1 mean free time. The dashed lines are Gaussian distributions with the same mean and same 
variance. These distributions have been obtained with ~4x 10 10 independent values of W(t). 



V. CONCLUSION 



In this paper we have given the details of the derivation of an equation for the cumulants generating function //(A) of 
the injected power in granular gases, in the case of a homogeneous random driving mechanism. This equation appears 
as a generalization of the Boltzmann equation for the velocity pdf of the gas. It has also the remarkable physical 
interpretation of being the kinetic equation of a gas of annihilating/cloning/colliding particles with an external energy 
source. This interpretation leads to our main result: the large deviation function ^^(w) of the injected power has no 
negative branch. We have also shown how to exploit a Sonine expansion to get better and better approximations of 
//(A) near A = 0. We finally obtained /i(A) in a closed analytical form by means of a generalized Gaussian hypothesis 
which is duly justified. This approximation becomes exact in the large dimension d — > oo limit. We have also 
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FIG. 14: p(w,t) (at t equal to 1 mean free time) divided by a Gaussian with same average and same variance for different 
values of a (at fixed constant temperature T g ) in the Driven Inelastic Maxwell Model. The light dashed lines represent a fit 
with a polynomial of fourth order. 



presented an argument that suggests the equivalence between the large deviation function of W(i) and that of the 
energy dissipated in collisions T>(t). All these results are consistent with the absence of negative events in Tr ao (w). 
This is at odds with previous studies that concluded that the GCFR could be satisfied for the integrated injected 
power in granular gases. Numerical simulations, on the other hand, show that the times at which a check of GCFR 
can be investigated are much smaller than the times required to reach the asymptotic large deviation scaling. At such 
small times the pdf of W(t) near W(t) =0 is almost Gaussian and this automatically leads to an artificial verification 
of the GCFR. On the other hand, the simulation of the inelastic Maxwell model instead of the inelastic Hard Spheres 
model made us able to gather much larger statistics so that we could see that even at finite times, for which irt(w) 
still displays negative events, the GC relation can be violated. 

We emphasize here that observing non-Gaussian features in the tail of the distribution of injecting power does not 
guarantee that the asymptotic scaling regime has been reached. Non Gaussianities may equally well arise from short 
time effects, and it turns out that the latter mechanism is at work already at very short times. Observing an apparent 
GCFR at such short times appears to be an artifact, and much larger times are required to sample correctly the 
asymptotic regime. In this respect, the time behavior of rescaled cumulants (as displayed in Fig. Ill|l is a valuable 
tool to decide if the "transient" regime is over. 

The analytical method put forward here to compute large temporal deviations is quite general and is particularly 
relevant in the context of the computation of probability distributions of global observables, which are useful, as 
underlined in the introduction, for the determination of universal features of non-equilibrium systems. For example, 
it applies to ballistically controlled processes, be they at equilibrium or not. Future work along these lines also includes 
the consideration of out of equilibrium granular gases without a steady state, such as freely cooling systems. 
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FIG. 15: (Color online). Finite time test of Gallavotti-Cohen relation for the injected power (with t equal to 1 mean free 
time), i.e. 7Tt(to) — 7ft (— w) vs. w, in a numerical simulation of the Driven Inelastic Maxwell Model with N = 50, and different 
values of a (the driving amplitude P has been rescaled in order to fix the granular temperature T g ). The dashed curve is a 
straight line with slope f3 = 1/T g . The dotted curve is a straight line obtained fitting the a = 0.1 data points until w = 45, 
useful as a guide for the eye. The thin (red) solid curve is a fit with a cubic (0.28ii) + 5.6 • W~ 4 w 2 — 1.1 ■ 10~ 5 w 3 ). 
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